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Abstract. Radiative acceleration of newly-formed dust grains and transfer of momentum from the dust to the gas 
plays an important role for driving winds of AGB stars. Therefore a detailed description of the interaction of gas 
and dust is a prerequisite for realistic models of such winds. In this paper we present the method and first results 
of a three-component time-dependent model of dust-driven AGB star winds. With the model we plan to study the 
role and effects of the gas-dust interaction on the mass loss and wind formation. The wind model includes separate 
conservation laws for each of the three components of gas, dust and the radiation field and is developed from an 
existing model which assumes position coupling between the gas and the dust. As a new feature we introduce a 
separate equation of motion for the dust component in order to fully separate the dust phase from the gas phase. 
The transfer of mass, energy and momentum between the phases is treated by interaction terms. We also carry 
out a detailed study of the physical form and influence of the momentum transfer term (the drag force) and three 
approximations to it. In the present study we are interested mainly in the effect of the new treatment of the dust 
velocity on dust-induced instabilities in the wind. As we want to study the consequences of the additional freedom 
of the dust velocity on the model we calculate winds both with and without the separate dust equation of motion. 
The wind models are calculated for several sets of stellar parameters. We find that there is a higher threshold in 
the carbon/oxygen abundance ratio at which winds form in the new model. The winds of the new models, which 
include drift, differ from the previously stationary winds, and the winds with the lowest mass loss rates no longer 
form. 
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1. Introduction 

The extended atmospheres of AGB stars are sites of large 
local and global variations in physical quantities which 
reflect the variability of the stars. It is here that dust 
grains icondense from the gas phase. The opaque dust is 
pushed out by the radiation pressure from the luminous 
star. The effects of interactions involving dust and shock 
waves caused by stellar pulsations can at the right condi¬ 
tions lead to the formation of a massive stellar wind. This 
process is critical for the evolution of AGB stars. The stel¬ 
lar wind does not only limit their lifetime but also enriches 
the surroundings with processed matter. First a circum- 
stellar envelope is formed and later the products are mixed 
with the interstellar medium. Observed long-time varia¬ 
tions in mass loss rates indicate that the properties of the 
stellar winds change as the stars evolve. The AGB phase 
ends when the star has almost completely lost its envelope 
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and soon afterwards appears as a white dwarf surrounded 
by a planetary nebula. 

Episodic mass loss variations of AGB stars on long 
time scales, i.e. 10 4 -10 5 years, have been observed in the 
form of detached CO shells for a long time (e.g. Olofssou 


et al. 1996 , 2000| ). Steffen fe Schbnberner ( 2000 ) argue 


that thermal pulses likely are responsible for the origin 
of the detached CO shells. Their circumstcllar envelope 
model includes separate equations of motion for gas and 
dust that are coupled to radiative transfer. In this context 
we also mention that there are variations on shorter time 
scales, on the order of 10 2 -10 3 years, that are believed to 
be associated with the duration of a helium shell flash at 
the beginning of a thermal pulse cycle. 

Mass loss variations on time scales of 10 2 -10 3 years 
that unlikely can be explained by thermal pulses have 
more recently been observed in the form of concentric 
arcs (concentric shells; e.g. Mauron & Huggins 1999 


mx and references therein). |5imis et al| ( [2001| . hence- 
forth SID01) draw the conclusion that a two-fluid gas-dust 
interaction produces mass loss variations on a time scale 
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of about 10 2 - 10 3 years, that is seen to agree with obser¬ 
vations of the dust-enshrouded star IRC +10216. Time 
dependent dust formation is included in their model. 

From a different perspective another physical mech¬ 
anism is proposed to play the key role in wind models 
of Late-AGB and Post-AGB objects. Boker ( [2002 ) argues 
that the concentric arcs (M-arcs) observed around these 
objects are unlikely to originate in the wind acceleration 
zone through the interaction of gas and dust. Instead, they 
could be the result of an (ad hoc) solar-like magnetic ac¬ 
tivity cycle in the star ([Boker| 200C|). Garcia-Segura et al. 


( | 2001 | ) also find, without including the dynamic effects 
of the dust component, that a solar-like magnetic cycle 
without mass loss variations reproduces many properties 
of observed concentric arcs. 

These studies motivate a closer investigation of the ef¬ 
fects of the dust-gas coupling on AGB wind structures. 
Not only is a closer study of the origin of the shells inter¬ 
esting. From a more fundamental point of view the gas- 
dust coupling is essential for the radiative driving of a 
stellar wind. It is important to study the limits of this 
coupling. Moreover, a model with improved physical ca¬ 
pabilities will provide the grounds for both qualitatively 
and quantitatively better estimates of mass loss rates and 
spectral energy distributions. The conditions for stellar 
dust formation can also be better understood. 

To describe the wind correctly the models have to in¬ 
clude a sufficient treatment of all three interacting com¬ 
ponents: gas, dust and the radiation field. Existing AGB 
wind models are either stationary or time-dependent. The 
models based on a stationary formulation do not admit 
flow variations with time. To their disadvantage few stel¬ 
lar parameter configurations have been shown to support 
stationary outflows (winds). On the other hand, time- 
dependent models tend to have (over-) simplified descrip¬ 
tions of radiative transfer (e.g. a semi-analytical treat¬ 
ment, or inadequate (often gray) opacities). 

Another model subdivision can be made regarding the 
degree of coupling between the gas and the dust compo¬ 
nents. In models assuming complete momentum coupling 
all radiative momentum gained by the dust immediately 
is transferred to the gas. Position coupled models, in ad¬ 
dition to complete momentum coupling, assume that the 
dust is mechanically bound to the gas phase, i.e. that it 
moves at the same velocity. 

The latest group of models in the literature however do 
not put any of the mentioned restrictions on the dust ve¬ 
locity. The degree of coupling inevitably affects the phys¬ 
ical distribution of both the gas and the dust in the en¬ 
velope. However, without detailed modeling it is not clear 
quantitatively how large the effect due to the coupling will 
be. The most recent works concerning the influence of the 
treatment of the gas and dust phases have been carried 
out by [Liberatore et al. (2001), SID01 and Steffen et al. 


( 1998| ). An overview of the handling of the gas-dust inter¬ 
action in earlier AGB wind models is presented by SID01. 

In this work we describe the method of building a 
three-component AGB star wind model. This model uses, 


to our knowledge, the most complete time-dependent de¬ 
scription of all three components, gas, dust and radiation, 
and it does not assume complete momentum coupling or 
position coupling. The physics and basic equations are pre¬ 
sented in Sect. |j. The numerical method and a discussion 
of different numerical approximations of the phase inter¬ 
action terms are given in Sect. |s]. This work is a study of 
dust-induced dynamic instabilities in the atmosphere. We 
study the dynamics of three-component wind models and 
compare the results with corresponding models where po¬ 
sition coupling, and hence complete momentum coupling, 
are assumed. The emphasis of the study is put on the 
effects of detailed momentum transfer. Section |l| contains 
the results and a discussion while the conclusions are given 
in Sect. |J. 


2. Physics of the wind model 


2.1. Original model characteristics 

The present work is based on an AGB wind model 


of Hdfner et al. ( 1995 , henceforth HFD95). We shall first 
summarize the properties of that model in this subsec¬ 
tion and then discuss the modifications in Sect. 2.2. The 


two-component (radiation hydro dynamic) RHD-sys tem of 
the stellar model as described by Feuchtinger et al. ( 1993 ) 
was combined with dust as a third component by HFD95, 
defining the RHDD-system. Some assumptions made in 
the RHDD-system concerning the dust and the gas are 
important in the current work and require some explana¬ 
tion. 

The matter in the wind is present in either of the two 
phases of dust or gas. The gas represents the overwhelm¬ 
ingly largest part of the matter (in some parts of the at¬ 
mosphere there is even no dust at all). The hydrodynamic 
equations describing the gas phase are the equation of con¬ 
tinuity, the equation of motion and the equation of (inter¬ 
nal) energy. A perfect gas law is adopted for the equation 
of state; the ratio of the specific heats 7 = 5/3, and the 
mean molecular weight g = 1.26. 

Only those particles in the gas that are part of the 
dust chemistry can move between the two phases. The 
dust phase is assumed to be composed of spherical dust 
grains in the form of amorphous carbon. The dust equa¬ 
tions that correspond to the equation of continuity for the 
gas phase (and describe the formation and destruction of 
dust grains) are the four moment equations for the mo¬ 
ments K 0 -K 3 of the grain size distribution function (Gail 
|fe Sedlmayi 1988 ; Gauger et al. 199C ). Dust formation 
is hereby treated self-consistently including the processes 
of nucleation, growth, evaporation and chemical sputter¬ 
ing (by gas particles) in a collision-less dust medium. The 
moments are related to the average of powers of the dust 
grain radius and allow the calculation of average proper¬ 


ties of the dust grains (Gail et al. 1984) such as: the total 
number density of dust grains rid = Kq', the mean grain 
radius (r^) = roKi/Ko (where ro is the monomer radius); 
the mean grain surface area (A) = Kq] the total 
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number density of monomers condensed into grains K 3 
gives the grain size ( N) = K 3 /K 0 ] the dust mass den¬ 
sity is pa = TO 1 -K 3 (mi is the dust grain monomer mass). 
The number densities of the gas-phase molecules that are 
involved in the grain formation are calculated in an equi¬ 
librium chemistry of H, H 2 , CO, C, C 2 , C 2 H and C 2 H 2 . 
The last four species contribute to the grain formation 
processes. All abundances are solar except for the carbon 
abundance which is specified through the carbon to oxy¬ 
gen ratio (ec/eo)- 

The model assumes complete momentum coupling; all 
momentum gained by the dust from the radiation held is 
immediately transferred to the gas. Assuming a very ef¬ 
ficient mechanical (position) coupling of the dust to the 
gas ( Dominik et al. 1989| ), the two phases are defined to 
move at the same velocity. The resulting equation of mo¬ 
tion is, 


d 

gj(pu) + V • {puu) = 

/pressure,g T /grav,g T /rad,g T /rad,d — 

Gm r 47 t 47 t . . 

-VP- 5 — p-\ - n„pH-\ -Kd pH (1) 

C ° C 

where /pressure,g is the force due to the gradient of the 
gas pressure; / gr av ,g is the gravitational force acting on 
the gas; / ra d,g is the radiative pressure force acting on the 
gas; / r ad,d is the radiative pressure force acting on the 
dust. Furthermore p is the gas density; u the gas (and 
dust) velocity; K g the (gray) gas opacity; Kd the (gray) 
dust opacity; P the gas pressure; H the first moment of 
the radiation field. Note that since Pd <C P ~ and therefore 
/grav,d <C /grav.g ” it is assumed that the only important 
term related to the dust in this equation is the radiative 
pressure acting on the dust. 

The dust (internal) energy equation is replaced by a 
radiative equilibrium relation since the dust is effectively 
thermally coupled to the radiation held. By assuming ra¬ 
diative equilibrium and LTE the dust temperature Xd is 
equal to the radiation temperature X ra d in the gray case. 
The radiation temperature is in turn defined by the radia¬ 
tive energy density J, 


T _ m 4 
^ ^ 1 rad 

7r 


( 2 ) 


where ctb is the Stefan-Boltzmann constant, J is the ze¬ 
roth moment of the radiation held. In addition the dust 
holds a negligible thermal energy compared to the radia- 
tive energy and the g as internal energy (cf. section 2 . 2 f in 
Hofner & Dorh 1992| , and references therein). 

The radiation held is described by the frequency- 
integrated zeroth and first moments of the radiation inten¬ 
sity. These moments represent the radiative energy density 
and radiative energy flux respectively. The corresponding 
moment equations of the radiative transfer equation are 
solved together with the hydrodynamic equations for the 
gas and the dust. At each time-step the equation of ra¬ 
diative transfer is solved for a given structure using the 


method of characteristics (e.g. Yorke 1980; Balluch 


This solution makes it possible to calculate the Eddington 
factor and other quantities that are necessary to close the 
moment equations. The gas and dust opaciti es are both as - 
surned to be gray. The gas opacity is set to (Bowen 1988), 


K g = 2 x 1CT 4 [cm 2 g" 1 ] 


( 3 ) 


while the dust opacity is (cf. Fleischer et al. 1992), 


Kd — Q ext. ^ 3 ; and Q ext 


Qext 

(m) 


( 4 ) 


where Q ext is the grain extinction efficiency. Previously, 
the Rosseland mean of Q axt was assumed as 
<3e Xt [cm -1 ] = 5.9 Td (used in e.g. HFD95; Xd is 
given in K). In later works (including this) it has 
been replaced with a better fit to the opacity data, 
<5e Xt [cm _1 ] = 4.4Xd (Winters 1994, priv. comm.; cf. 
Hofner & Dorfi 1997). The dependence of Q' ext on Xd 
follows as a consequence of taking the Rosseland mean 
of the frequency-dependent extinction coefficient. The 
consistent treatment of the radiation held is a strength 
of the RHDD-system model compared to models using 
semi-analytical approximations. A more recent formula¬ 
tion of th e RHDD-sys t em also includes non -gray radiative 
transfer ( Hofner 1991; Hofner et al. 2002a ,[]). 


One physical limitation of the RHDD-system descrip¬ 
tion in previous papers is the assumption of position cou¬ 
pling. By this assumption the effects of two drifting phases 
are totally disregarded, and can therefore not be properly 
considered. A separation of the two phases does not only 
require an additional equation of motion for the dust com¬ 
ponent but also several phase interaction terms. 

The naming convention in the literature of AGB wind 
models varies depending on the description of the radia¬ 
tion held or the presence of a freely moving dust compo¬ 
nent. Models in which the focus is on dynamics are often 
called “n-huid” models where the “n” denotes the num¬ 
ber of separate equations of motion for different material 
components. We prefer to label the models according to 
for how many physical components the conservation laws 
are solved - counting both material phases and the ra¬ 
diation held emphasizing that we neither use simple 
equilibrium assumptions, nor prescribe values for a par¬ 
ticular component which appear in interaction terms of 
other components. In this sense, we refer to our models as 
three-component models (gas, dust and radiation held), 
not as single- or two-fluid models. Furthermore, in our 
notation the RHDD-system models are position coupled 
(PC) three-component models 1 ] as opposed to the three- 
component drift models that are described in the following 
subsection. 


1 In spite of the PC in these models the dust is still treated 
as a separate component, since the time-dependent formation, 
growth and evaporation of grains are described by the moment 
equations. 
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2.2. The dust equation of motion 


Most studies of dust-driven stellar winds have evolved 
around the assumption of complete momentum coupling 
in stationary winds. They also often contain a simplified 
description of either the radiation field or the dust com¬ 
ponent (e.g. instantaneous dust formation and a constant 
grain size). The two latest works on drift in stellar winds of 
AGB stars are those by Liberatore et alj ( [2001 ) and SID01 
(we refer to SID01 for an overview of previous studies con¬ 
cerning the effects of drifting phases in cool stellar winds). 
The former have carried out a study of different degrees 
of dust-gas coupling in stationary winds of late-type stars 
using a frequency dependent dust opacity. However, the 
grain radius in the dust component is assumed constant. 
SID01 have carried out explicit time-dependent hydrody- 
namical modeling in a two-fluid medium using a given 
temperature structure, but treat dust formation in detail. 

In this study we relax the assumptions of position cou¬ 
pling as well as complete momentum coupling made in the 


previous modeling (see Sect. 2.1). We do this by adding an 


equation of motion for the dust component to the RHDD- 
system and modifying the equation of motion for the gas 
accordingly. Phase interaction terms are included to con¬ 
serve the physical quantities (see the following subsec¬ 
tion). The resulting system is henceforth referred to as 
the RHD3-system, where the third “D” stands for drift; 
the models are accordingly named drift models. 

The dust component consists of dust grains of differ¬ 
ent sizes. In principle each group of dust grains of a cer¬ 
tain radius can be ascribed a separate equation of motion. 
Coupling terms between the equations for dust grains of 
different sizes may be neglected on the assumption that 
grain-grain collisions are far less frequent than grain-gas 
collisions (the dust is assumed to be pressure-less). We can 
then take a mean of the velocity equations of individual 
sizes to get the size-averaged dust equation of motion. This 
mean equation formally looks like an equation of motion 
for one grain size. 

In the new system the gas equation of motion, Eq. (I:), 
is exchanged with, 


d 

— (pw) + V • {puu) = 

CtTYI 4l7T 

- VP- -jTp + —Kg pH + / drag - ScondW 1 (5) 

where the last two terms represent the momentum trans¬ 
fer by gas-dust collisions and the momentum change 
when dust condenses from the gas phase respectively (see 
Sect. 2.3.1 ). The pressure-free dust equation of motion in 
turn is, 


d_ 

dt 


{p d v) + V • {p d v v) = 


2 Pd + S',] p 11 /drag T *^cond a (b) 


v is now the “mean” dust velocity and pd is the dust den¬ 
sity. 


By the same reasoning as in the RHDD-system we do 
not include the dust equation of internal energy but as¬ 
sume that the grain temperature is determined by radia¬ 
tive equilibrium. 


The dust formation processes are affected in several 


ways by drift (e.g. I 

'Huger & Sedlmayr 

1997, 

henceforth 

KS97; 

Dominik et a; 

. 1989; 

Draine & Salpetei 

1979 

). We 


do not include these modifications in the models presented 
in this paper, but plan to do it in the future. 


2.3. Phase interaction terms 

In the gas-dust phase interaction we must consider each 
of the three transferred physical quantities of mass, mo¬ 
mentum and internal energy. There are two types of mo¬ 
mentum and energy exchange between the gas and dust 
phases. On the one hand they are transferred with the 
mass that switches phase, and on the other hand they are 
transferred in the collisional interaction. 

2.3.1. Mass transfer interaction terms 

Mass is transferred between the phases when dust grains 
form or grow or alternatively when they evaporate. The 
rate at which material condenses iS con d/mi corresponds to 
the r.h.s. combination of source terms in the Kj, equation 
(Eq. 7 in HFD95). The rate <S CO nd ' s also a sink term in 
the gas equation of continuity. The rate of the momentum 
transfer at mass transfer is represented by the last term 
in Eqs. (|| & ^). We set the velocity of the formed (or 
evaporated) dust u 1 equal to the gas velocity u. The rate 
at which internal energy is transferred, and work is done 
by removal of mass from the gas phase, S con dh n , is a sink 
term in the gas internal energy equation. Here h n is the 
specific enthalpy, 

h n = e n + P n /p n . (7) 

The superscript n indicates that, on addition of mass to 
the gas, these quantities might be in non-equilibrium with 
the gas. 

The effects of the mass transfer on the gas are very 
small compared to the effects on the dust. Nevertheless 
we include all terms above in both the dust and the gas 
equations for completeness, with the exception of the in¬ 
ternal energy transfer term (<S C ond^ n )- The internal energy 
transfer to (and from) the gas phase is ignored on the basis 
that the internal energy transfer with the radiative field 
always will be (several) orders of magnitude larger. 

In contrast to the terms discussed next all interaction 
terms mentioned so far are independent of the presence of 
the separate dust equation of motion 

2.3.2. Collisional interaction terms 

The momentum transfer term, the drag force, is the most 
important interaction term in dust-driven winds. Most of 
the radiative pressure acts on the opaque dust grains ac- 
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celerating them outwards. The gas is dragged along by the 
accelerating dust and forms a stellar wind. 

The drag force is derived from the local physical con¬ 
ditions. In our case these can be characterized as follows. 
The gas particle velocities have a Maxwellian distribu¬ 
tion (the gas is described with the continuum approxima¬ 
tion) . The dust grains in the dust medium are in the free- 
molecular regime compared to the gas. In App. [A| we dis¬ 
cuss the validity of this and the previous assumption in the 
current context. The dust grains are assumed to be spher¬ 
ical. Collisions between gas particles and the dust grain 
surface can be either specular or diffusive, depending on 
how the normal and tangential momentum is distributed 
in the collision. In a specular collision of the incident par¬ 
ticle the normal component of the velocity (in the frame 
of reference of the dust particle) is reversed on reflection, 
while the tangential components are unaffected. In a dif¬ 
fusive collision the particle is first acommodated on the 
surface, then it is thermalized, and finally it is emitted in 
a random direction. The drag force is derived by integrat¬ 
ing the pressure over the surface of the dust grain (cf. e.g. 
Hayes fe Probstein 1959 ; Schaai 1963| ). The general form 
of the drag force is 


It is not yet established how collisions are distributed be- 


,■ _ e D 

J drag — (JpTl d 


C D 


a = 7T(r d ) 2 = /Kq . 


(9) 


Cd = Co,diff + Cd, 


2 7T2 (1 — e) 


common — 


+ 


+ 


4'S'd 


■4SS-1 


2 ^ 


erf(Sb) 


2<S'q + 1 , c 2, 

i c3 exp(-S D ) 
71-2 5 d 


( 10 ) 


c , 2 k B T g 

od =-, where u mp = W- 

Vmp V ^ TOH 


( 11 ) 


tween specular and diffusive (see however e.g. Kruger et al 


1994, henceforth KGS94), and we leave the option to use 


( 8 ) 


where a denotes the gas-dust geometrical cross section; v B 
is the drift velocity (v — u)\ Cb is the drag coefficient. The 
factor of 2 comes from the definition of the drag coefficient. 
For a we use, 


Because of its concisene ss we prefe r the form of the drag 
coefficient presented by Bird (|1994, henceforth B94), 


where erf is the error function; T g is the gas kinetic tem¬ 
perature; the fraction of specular collisions is defined by £ 
(see below); Sb is the speed ratio, 


a combination of them. 

In the collisional interaction the internal energy of the 
gas is modified by on the one hand inelastic (diffusive) 
dust-gas collisions (cf. the discussion on ‘g acc ’ in KGS94). 
On the other hand kinetic energy is converted into inter¬ 
nal (thermal) energy as gas particles which preferentially 
come from one direction (the drift) are reflected in ran¬ 
dom directions when hitting a dust particle. The energy 
transferred to the gas in this process, in addition to the 
work done by the drag force on the gas, corresponds (in 
the case of only specular collisions) to the relative speed 
of the gas and dust particles times the drag force (cf. the 
treatment of ‘<?f r i c ’ in KGS94). We ignore the effects of 
both these terms ( q acc and gf r i c ) in the models presented 
in this article, on the assumption that the effects on the 
gas internal energy are minute when compared to the ra¬ 
diative energy exchange with the gas. Preliminary results 
of our time-dependent models including the second heat¬ 
ing term (<7f r i c ) show that its influence on the structure is 
negligible. 

In connection to this discussion we want to point out 
that the same term was found by KGS94 to play a sig¬ 
nificant role in the gas energy balance in their stationary 
models. It is, however, difficult to compare their results 
with ours because the models differ in several respects. 
The most important differences are: different physical as¬ 
sumptions on the radiative transfer; different stellar pa¬ 
rameters of the models; they use constant sized dust grains 
(i.e. no time-dependent formation). We plan to discuss the 
importance of this term in a forthcoming article. 

2.4. Complete momentum coupling 

The validity of the assumption of complete momentum 
coupling (CMC) in the cool stellar dust-driven wind has 
been subject of extensive discussions starting with Gilman 
(1972). For CMC flows it is assumed that all radiation 
momentum is immediately transferred to the gas, and the 
inertia of the dust phase is neglected. Hence, 


./'drag ,/rad,d T ,/grav.d 


( 12 ) 


v mp is the most probable thermal speed of the Maxwellian 
velocity distribution. 

The second term on the r.h.s. in Eq. (0), Cb, common 
(the common term), is sufficient when all collisions are as¬ 
sumed to be specular. The first term, Cb,diff i is in addition 
required in diffusive collisions. £ = 1 corresponds to fully 
specular collisions and £ = 0 to fully diffusive collisions. 
A combination of collisions is achieved by using values on 
£ in between (however as pointed out by B94 this combi¬ 
nation is only an approximation to a real scattering law). 


where / gr av ,d is the gravitational force acting on the dust. 
With only specular collisions - and using one of the ap¬ 
proximative for ms of the drag coefficient Cd for /drag pre¬ 
sented in Sect. 3.2 - this expression can be inverted with 


respect to the (equilibrium) drift velocity Td (see e.g. sec¬ 
tion 2.2.3 KS97). With this relation there is no need to 
solve the dust equation of motion, and computational time 
is saved when solving the system of equations without it. 
However when we consider diffusive collisions an inver¬ 
sion is not possible anymore and Eq. © must be solved 
numerically. In that case there is no computational mo¬ 
tivation to use the assumption of CMC, and one may as 
well solve the dust equation of motion. We do not assume 




























6 


C. Sandin and S. Hofner: Three-component modeling of C-rich AGB star winds 


CMC in the RHD3-system but we give a qualitiative cri¬ 
terion on how close our models are to CMC in Sect. 4.3. 
Equilibrium drift expressions are used by e.g. KS97 and 
in some of the calculations presented in SID01. 


3. Numerical method 

Before we look at the numerical differences of the RHD3- 
system compared to the RHDD-system we summarize the 
main features of the RHDD-system. 


3.1. Features of the RHDD-system 


A detailed description of the numerical method of the 
RHDD-system is found in Dorfi & Hofner ( |1991 ). The 
RHD-system, without dust, was described in Dorfi & 
Feuchtinger ( 1995| , |1991| ). 

The gridpoints are distributed with an adaptive 


grid (Dorfi & Drury 1987) in which a grid equation re¬ 
solves gradients of selected quantities. Currently the grid 
resolution function is determined by the thermal (inter¬ 
nal) energy and the gas density. The temporal smoothing 
factor is set to r g = 10 2 s, which is orders of magnitude 
smaller than, e.g., dynamical or dust time scales in the 
problem, meaning that the grid can freely adapt to phys¬ 
ical features; the spatial smoothing factor is set to a = 2. 


Artificial tensor viscosity ( Tscharnuter fe Winkler 


1979|) is used in regions subject to inhomologous contrac¬ 


tion. The term is added as a source term in the gas equa¬ 
tion of motion and the (internal) energy equation. The 
shock front is thereby widened to the relative characteris¬ 
tic length scale l, 


3.2. Numerical issues in the RHD3-system 


In this section we address several numerical issues associ¬ 
ated with the dust equation of motion. The RHD3-system 
now consists of thirteen equations. 

The grid equation can be adjusted to resolve both the 
gas and the dust components by including corresponding 
dust qua ntities in the grid reso lution function (as sug¬ 
gested by Dorfi fe Gautschv 199(]| ). If this is done, however, 
the number of gridpoints should be increased to resolve 
both components. An increased number of gridpoints has 
the disadvantage that the fraction of mass contained in 
some grid cells may become smaller than the numerical ac¬ 
curacy of the scheme, thereby introducing new problems. 
Hence we leave the grid resolution function unchanged 
(compared to HFD95). Presently we use 500 gridpoints in 
all calculations, and a first order donor cell advection is 
adopted in all drift models. 

An artificial viscosity term analogous to the one in the 
gas equation of motion is added to the dust equation of 
motion. Like the gas shocks, strong dust velocity gradients 
are widened by a characteristic dust front length scale (see 
Sect. 4.2 for a description of the term dust front), i.e. 


Id = rfd 


(14) 


We set both / and /d equal to 3.5 10 -3 in all our calcula¬ 
tions. 

In drift models dust density gradients at shock 
fronts may become extremely steep, despite the widening 
achieved with the artificial viscosity. To smear out these 
gradients we add artificial diffusion in these models in the 
form presented by WN86. This term is added as a source 
term in all four dust moment equations, 


D k , = V • (skVIQ 


(15) 


l = rf 


(13) 


where r is the local scale length, i.e. the radial distance 
from the center, and / is a constant that defines the width 
of the shock front as a fraction of r. 


In addition to the five RHD-equations, the four dust 
moment equations and the grid equation there are two 
more equations; an equation of the integrated mass and an 
equation keeping track of the condensible amount of car¬ 
bon. Thus totally there are twelve non-linear equations, 
out of which ten are partial differential equations (PDEs). 
All equations are discretized in the volume-integrated con¬ 
servation form on a staggered mesh (Winkler & Norman 
henceforth WN86). The spatial discretization of the 


advection term can be chosen to be either first order 
(donor cell), or second order (monotonic advection 


vail 


Leer 1977). The same order of precision is used in all 
PDEs. 


The full RHDD-system of twelve equations is solved 
implicitly using a Newton Raphson algorithm where the 
Jacobian of the system is inverted by the Henyey method. 


where Ki denote the dust moments Kq-K^. The transport 
coefficient c is in general defined as, 

S' = l 2 /r (16) 

where l is defined as above. We define the characteristic 
shock propagation time r, which is the time needed to 
cross a region of relative width Ax as, r = Aai/lwj, where 
w is the shock velocity. With the dust velocity s as a 
measure of the shock velocity and the local scale length r 
representing the relative width Air we have, 

Sk = flr\v \. (17) 

A physically motivated explanation for the use of arti¬ 
ficial mass diffusion in our drift models can be found in the 
experience from other areas of numerical hydrodynamics. 
As WN86 (and references therein) discuss, spurious results 
occur when strong shocks interact with walls, contact dis¬ 
continuities and other strong shocks. Normally, there is 
enough numerical diffusion implicit in the advection of 
the numerical scheme to prevent these features from ap¬ 
pearing, but if and when there is not they may appear as 
strange spikes in the solution. 
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In our drift models features (spikes) that can be 
attributed to the interaction between steep gas shocks 
and dust fronts appear first in the dust velocity, and 
subsequently in the dust density (i.e. the dust moments). 
Furthermore the time steps are decreased by the large 
variations of these quantities. Artificial diffusion prevents 
most of these features from appearing, but does not 
manage to remove all of them (see Sect. [h^). The model 
evolution and the wind characteristics are, however, not 
affected by their presence since the dust density in the 
region of the feature always is very small. While the de¬ 
scribed features are interpreted by the artificial viscosity 
term as contracting regions on the outwards facing side 
(i.e. away from the star), they are not interpreted as such 
on the inside. Consequently, these features can always be 
identified by the “discontinuous” jump on the inwards 
facing side. When we make estimates of the maximum 
dust speed, dust density variations and related quantities 
we must first separate these features from the rest of the 
structure: - 


f — 
Jcond — 


K, 


K.i + n c 


<e = / ( 


min 

cond 


(18) 


12 x 10" 4 ~ 10“ 3 . 


(19) 
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Fig. 1. Relative error of the approximative drag coeffi¬ 
cients with respect to the expression Cb (Eq. 0 ). The 
limits approximation (Eq. pis| ) is represented by the three 
lines grouped at Crt A ; the approximation of Berruyer & 
Frischl ( |1983| ) (Eq. g) is represented by the three lines at 


The dust equation of motion (Eq. |(j) is not numerically 
well defined in regions with very small amounts of dust. 
On the other hand, very small amounts of dust are not 
likely to affect the structure of the stellar wind. The dust 
moments become irrelevant for low degrees of condensa¬ 
tion (/cond) that satisfy, 


where p / ot is the total density of condensible matter (both 
that present in the gas and the dust phases); n c is the 
total number density of condensible material in the gas 
phase; e is the numerical limit of an insignificant amount of 
dust. From our numerical experience we have found e to be 
about 30 / . The total abundance of atoms of condensible 


Cp 4 ; the single line C^ v represents the high-velocity ap¬ 
proximation (Eq. ^4|). The discontinuity in Cq V at (1,1) is 
due to the discontinuous plot ordinate. Solid lines are com¬ 
puted assuming fully specular collisions (e = 1). The dot¬ 
ted lines assume fully diffusive collisions (e = 0). For the 
dash-dotted lines 50% of the collisions (e = 0.5) are spec¬ 
ular and the rest diffusive. The temperature ratio adopted 
in the diffusive term Cb.diff is 1.0. The vertical bars show 
the range of typical values of the speed ratio Sb in our 
wind models (see Eq. p5|). 


The part of the analytical drag coefficient common to 
both specular and diffusive collisions, Cb, common (Eq. pT 


materi al io about 10~ 4 of the total number of gao particles 
(this is the carbon that is not bound to CO). Thus for 
carbon, 


is numerically difficult for small speed ratios Sb- Bainei 
et al. (1965|), Berruyer fe Frisch (1983) and others have 


We find the numerically limiting quantity of “no dust” to 
be, 


derived approximations to this expression. Here we com¬ 
pare three approximations of the drag coefficient with the 
expression in Eq. [HJ All three approximations are writ¬ 
ten below in the same form as Cb- We have added the 
diffusive collisions term to the first two approximations. 
[Berruyer fe Frisch] derived the following approximation of 
Cb, common (see also, Liberatore et al. 2001; Mastrodemo; 

Tf§), 


et al 


Pd 

P 


Pc 


•/< 


min 

cond 


10 “ 3 • 10“ 7 = 10 ~ 10 . 


( 20 ) 


no dust 


c£ F = 


In regions where there is dust and the dust/gas ratio even¬ 
tually becomes lower than the limit given above we ex¬ 
change the dust equation of motion, Eq. (|(ij), with the re¬ 
lation, 


-m- 


2 7T2 (1 — e) 

3 

1 + s^- 


C 2 


sgn(Sb) 


(0.17S& + 0.5|S D | + 1) 3 . 


( 22 ) 


In the limits approximation the separate solutions of high 
and low Sp are combined (see e.g. KG S94, MacGregor & 
Stencel 1992) ; Draine & Salpeter 1979), 


For reasons of stability we use this relation also in dust 
forming regions where the dust/gas ratio is smaller than 
10 times the value in Eq. (p0|). 
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The influence of the thermal conditions is ignored in the 
high-velocity approximation (see e.g. SID01), 


= 2 . 


(24) 


Figure [I] shows the relative errors of the approximative 
drag coefficients presented above compared to Cd- In the 
figure we see that C® F deviates the least from Cb- The 
relative error of Cq A is never larger than 1.5% and it is 
easier to code than C§ F and therefore we use C^ A in our 
calculations. 

In their work on a two-component wind model SID01 
used the high-velocity approximation of the drag coeffi¬ 
cient in the momentum transfer. This choice was moti¬ 
vated by the fact that the momentum transfer term makes 
the system of equations too stiff to be solved explicitly, 
thereby requiring very small time steps. To make an an¬ 
alytical time-integration of the drag force possible (and 
thereby allow larger time steps) t hey chose the simpler 
form achieved with Cq V . In Sect. T3 we study the re¬ 
sponse of the RHD3-system to the two different drag co¬ 
efficients Cq V and Cq A . 

The advantage of an implicit scheme is that the 
time step is not limited by the very restrictive Courant- 
Friedrichs-Lewy (CFL) condition. In our models it is in¬ 
stead limited by the temporal variations in the physical 
quantities. Currently local variations less than 20% be¬ 
tween two subsequent time steps are accepted or else the 
time step is reduced. Our experience shows that the dust 
velocity has the strongest variations of the quantities and 
therefore limits the time steps in drift models. The vari¬ 
ations of the dust velocity are the largest in gridpoints 
containing very little dust mass. A negative effect of a 
shorter time step is that variations on shorter time scales 
(noise) may occur in other quantities as well. Slightly more 
than ten times as many time-steps are needed to reach the 
same age in drift models compared to PC models. 

Note that the results presented here for the position 
coupled RHD3-system models do not exactly match the 
results of HFD95. The d evia tions are caused by a modified 
dust opacity (see Sect, [ 2 . l| ) and a few minor changes in 
how the dust moment equations are activated when the 
first dust grains form close to the photosphere. 


4. Results and discussion 

4.1. Modeling procedure 

The basic modeling procedure of all winds in this article 
is similar to that in HFD95. Note however that we use the 
RHD3-system code in all our calculations, including the 
PC models. 

The modeling procedure is as follows. All wind models 
are started from hydrostatic dust-free initial models where 
the outer boundary is located at about 2 R*. All five dust 
equations are switched on at the same time and dust starts 
to form whereby an outward motion of the dust and the 
gas is initiated. The expansion is followed by the grid to 
about 25 R* (typically 110 15 cm). At this radius the outer 


Table 1 . Stellar parameters of the models. The models in 
this article are named in correspondence with the models 
in HFD95, Table 2. 


model 

AT [Mq] 

L * [Lq] 

Teft [K] 

R* [Rq] 

A 

1 

to 4 

2600 

493 

B 

1 

to 4 

2500 

533 

C 

1 

to 4 

2400 

578 

D 

1 

3 10 4 

2600 

853 


boundary is fixed allowing outflow, and the model evolves 
for about 100 years (and more). The model calculation is 
stopped before a significant depletion of the mass inside 
the computational domain occurs. 

In the models presented here it is the dust that initi¬ 
ates and drives the stellar wind. The dust formation zone 
is always inside the model domain. An inflow of matter 
through the inner boundary which is located well below 
the photosphere (typically at about 0.9 R*) is not allowed. 
The fraction of the mass of the star contained in the model 
domain is about 15-30%, depending on the stellar param¬ 
eters. This large fraction is a consequence of usin g a rather 
low gas opacity for the pr esent models ( Sect. 
problem was discussed by Hofner et al. (|l998 , 


2 . 1 ), and the 
cf. Sect. 3.1 


and Fig. 1). Note that a too high density in the stellar 
photosphere does not necessarily translate into unrealistic 
conditions in the wind acceleration zone since the density 
structure in this region is strongly influenced by dynam¬ 
ical effects. A smaller mass fraction is achieved by using 
a more appropriate description of the gas opacity, which 
we plan to do in a forthcoming article. We do not make 
specific assumptions on the presence of drift between the 
gas and dust phases in different parts of the atmosphere. 

The model is determined by the three stellar parame¬ 
ters of stellar mass M, luminosity L and the effective tem¬ 
perature T e ff. In addition the carbon abundance is speci¬ 
fied through the carbon/oxygen ratio £c/eo• The model 
parameters in our study are given in Table [I]. The param¬ 
eters of the models we study in this section match those 
of HFD95, for clarity the models are named with the cor¬ 
responding letters. 

The discretization of the advection terms of the PC 
models are either first or second order in precision while 
those in the the drift models, always are of first order 
precision (see Sect. ||^). To estimate the possible influence 
of 1 st vs. 2 nd order numerics on the model we present 
results using both alternatives in the PC models. 


4.2. Three-component dust-driven wind models 

In this section we study the physical differences of drift 
models compared to PC models. In selecting the parame¬ 
ters for the models we wanted to cover both less massive 
and massive winds. What we describe is a dust-driven in¬ 
stability. Presently there are no pulsations since the use of 
a variable inner boundary condition that simulates pulsa- 
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tions could complicate the interpretation of other effects. 
Without an atmospheric levitation by (pulsational) shocks 
the models require higher ec/^o ratios to form winds. The 
results of the models are given in Table [|. Before we look 
at the results, we want to review the characteristics of the 
different winds that form. 

Stationary winds by definition do not show variations. 
In PC models it has been found that low mass loss winds, 
with low ec/eo ratios, can become stationary. With higher 
ec/eo ratios the dust formation process becomes unsta¬ 
ble and the models consequently time-dependent. In the 
back-warmed region behind a newly formed dust shell of 
such a model the temperature eventually becomes too high 
for further dust formation to take place, and the process is 
shut off. Radiation pressure meanwhile acts on the individ¬ 
ual grains to push the dust shell (and the gas) outwards. 
When the gas temperature behind the leaving dust shell 
has decreased sufficiently the dust formation process is re¬ 
activated and a new dust shell can form. This mechanism 


was called dust induced /{-mechanism by Fleischer et al. 


1995 and HFD95. 


The conditions of wind formation change in drift mod¬ 
els. The dust tends to accumulate in regions where the 
gas-dust interaction is strong. Strong interaction regions 
are typically represented by the dense gas behind shocks. 
When the interaction is not strong enough for the dust to 
drag the gas along, dust alone leaves the envelope and no 
significant wind forms. The dust component is a pressure¬ 
less medium that does not form ordinary shocks. However, 
dust tends to accumulate in high (gas) density regions 
which usually occur behind gas shocks. The result of the 
accumulation are strong gradients in the dust density and 
the dust velocity that coincide with the gas shocks. We 
will hereafter refer to these gradients as “dust fronts”. 

We discern between local and global differences in the 
wind. Global differences, such as that of the mean mass 
loss rates will be discussed after we have looked at local 
differences. Figure ^ shows the drift model C19dl and the 
PC model C19c2 at an evolved stage where the initial 
transient effects from the expansion phase have left the 
model domain. It is difficult to discuss periods of shock 
formation as C19dl develops into a wind showing irregular 
variations, the same is true for all drift models presented 
here. However, typical time scales are here on the order 
of years. At this point we want to stress that it is the 
general pattern we refer to in this and all following figures 
of the physical structure; it is difficult to compare specific 
features in winds showing irregular variations. The instant 
presented in Fig. has been selected because differences 
between the two models are seen clearly. The qualitative 
results that we discuss for this model are also valid for the 
other combinations of PC and drift models in Table ||. 

The dust density (Fig. |1) changes by about two or¬ 
ders of magnitude across dust fronts in C19dl, this is to be 
compared to the closer to one order of magnitude changes 
in C19c2. We can see a tendency towards narrower dust 
shells in C19dl compared to C19c2 in the dust-gas density 
ratio, Fig. and in the degree of condensation (Fig. |(h) 


we see that regions in front of the dust fronts are sharply 
depleted of dust. Without the position coupling constraint 
dust tends to leave the low density zones and is accumu¬ 
lated to the regions behind shocks. A consequence of the 
dust relocation is that dust is shifted to regions resolved 
by more gridpoints. The increased gridpoint density in the 
regions of dust fronts allows for a numerically more accu¬ 
rate calculation of the advection, which is crucial when 
using a first order numerical scheme. The correspondence 
of the first order drift models in terms of more numerous 
and sharper shocks is better when compared to the sec¬ 
ond order PC models than they are with the first order 
PC models. 

The dust velocity plot (Fig. |b) requires some extra 
explanation. As can be seen in the plot the dust has ac¬ 
quired e norm ous velocities at 15 and 25 i?*. As is discussed 
in Sect. T2 these features (also seen in Fig. i-j) appear 
in connection with steep dust fronts in drift models. Since 
they are always associated with very low dust densities 
they are not relevant for the model evolution and need 
not be taken into account in physical consideration. The 
cancellation of the high velocity vs. low density in the fea¬ 
tures is seen in the dust momentum density, Fig. |2|f. 

The drift velocity otherwise stays fairly low. The lo¬ 
cations of higher drift speeds are mostly the dust density 
troughs in the innermost region where the radiation pres¬ 
sure is the strongest. The maximum drift speed covering 
the calculated evolution of C19dl is about 30 kms -1 in 
inter-shock regions. Mostly, however, the upper drift speed 
limit is closer to 15 kms -1 (as is seen in the figure). Most 
of the dust is found in the regions behind shocks where the 
drift velocity is low by the increased number of collisions 
with the gas. We have found the lower limit of drift speeds 
to be about 0.1 kms -1 . The speed ratio most of the time 
stays within the limits, 


-l<log(|S D |)<l 


(25) 


The mass loss generated in C19dl is fairly large. In 
Table || we see that the massive winds of the C and D 
models are more easily formed in drift models than the less 
massive outflows. The evolution of the less massive winds 
formed in the B models are computationally more difficult. 
The models cover a shorter time interval and the resulting 
winds give less precise mean values. The parameter con¬ 
figurations that lead to the low mass loss time-dependent 
models and the stationary models in PC models (A with 
sc/^O < 2.5 and B with ec/eo < 2.0) do not result in a 
wind at all with drift included. The dust particles that ini¬ 
tially form in the supersaturated areas of the hydrostatic 
atmosphere do not get a chance to grow very large. In a re¬ 
peated sequence the dust first accelerates and brings the 
gas along. The removed position coupling together with 
the slow dust formation, however, prevent a further ex¬ 
pansion. The dust is not able to support the gas, which 
falls back while the dust leaves the atmosphere. The con¬ 
densible material in the hydrostatic initial model is not 
enough to drive a wind. A conclusion is that there is a 
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5 10 15 20 25 5 10 15 20 25 


Fig. 2. Spatial structures of the drift model C19dl (solid line) and the PC model C19c2 (dash-dot-dotted line): a gas 
velocity; b dust velocity; c gas density; d dust density; e gas momentum density; f dust momentum density; g dust-gas 
density ratio; h degree of condensation; i drift velocity; j dust-grid relative velocity. The features in the drift model 
associated with the inwards facing discontinuous peaks in the dust velocity and the dust-grid relative velocity (at 15 
and 251?*) are irrelevant for the evolution of the model and should not be taken into account in physical consideration 
(see Sect. and 4.2). 
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Table 2. Model quantities for the PC models (2nd- and lst-order precision in the advection terms; suffixes c2 and cl) 
and the drift models (suffix ell). The model names are constructed by combining the symbol in Table |l] representing 
the combination of stellar parameters with (10 x) the carbon/oxygen ratio. The models that we were not able to run 
for a longer time interval are marked with parentheses; the associated numbers are less reliable because the means are 
taken over small time intervals. The different types of winds are: s, stationary wind; i, irregular wind; q, periodic or 
quasi-periodic wind; —, no wind. 


model,ec/eo 

Type 

position cou 

2nd order (c2) 

(M) <Uoo) 

[ M ©yr _1 ] [kms -1 ] 

pled models 

1st order (cl) 
Type (M) 

[ M ©yi' _1 ] 

(Uoo ) 

kms -1 ] 

Type 

drift models 

1st order (dl) 

(M) 

( M 0 y f *1 

(Uoo ) 

kms -1 ] 

A23 

s 

2.0 10 -7 

11 

S 

2.3 10 -7 

ii 




A25 

i 

5.0 10 -6 

33 

s 

5.110 -7 

18 

— 



A27 

i 

4.4 10 -6 

35 

i 

6.4 10 -6 

35 

i 

(6.6 10 -6 ) 

(40) 

B20 

q 

9.7 10 -6 

16 

s 

2.0 10 -7 

16 

— 



B21 

q 

9.4 10 -6 

29 

i/s 

3.0 10 -6 

24 

i 

(i.i nr 5 ) 

(28) 

B22 

q 

1.0 10 -5 

34 

i 

9.6 10 -6 

27 

i 

(7.8 10 -6 ) 

(28) 

C18 

q 

1.4 10 -5 

24 

i 

1.3 10 -5 

21 

i 

1.5 10 -5 

21 

C19 

q 

1.5 10 -5 

26 

i 

1.4 10 -5 

22 

i 

1.4 10 -5 

24 

C20 

q 

1.5 10 -5 

30 

i 

1.5 10 -5 

24 

i 

1.3 10 -5 

25 

D16 

q 

5.3 10 -5 

28 

i 

5.6 10 -5 

23 

i 

4.3 10 -5 

27 


higher threshold in ecj^o when drift models form winds, 
compared to the PC models. 

From the numbers on the mean mass loss rates and 
mean terminal velocities in Table || we see that the drift 
models largely reproduce the same numbers as the corre¬ 
sponding PC models for those parameters that lead to the 
formation of a wind. The local differences previously dis¬ 
cussed therefore do not seem to affect the model globally 
in modifying the mean mass loss rate.^j 

The winds calculated with a first order numerical 
scheme do not appear periodic in any of the models we 
have studied. Periodicity puts very high demands on the 
numerical scheme. The first order scheme, maybe in com¬ 
bination with numerical noise allowed by the small time- 
steps (due to the gas shock-dust front interaction), pre¬ 
vents periodic models from forming by disturbing the 
models “randomly”. Compare the type of wind formed in 
first order precise PC models with those of second order 
in Table |2j. 

SID01 have found long-time modulations of several 
10 2 * yrs in their model that they attribute to the dust- 
gas interaction. Those variations (producing structures 
comparable to the observed concentric arcs) should how¬ 


2 The numbers in Table of the PC models do not exactly 
correspond to the numbers in HFD95, and we attribute these 

differences partly to the numerical differences of the codes we 

have used. In add ition, however, we use a different dust opac¬ 
ity (see Sect. 3.2) which probably accounts for most of the 
differences. 


ever not be confused with the variations in our model 
that are due to the dust induced ^-mechanism and hap¬ 
pen on much shorter timescales of typically a few years. 
The authors give the stellar parameters of their model 
as: M* = 0.7M 0 , L* = 2.41O 4 L 0 , T* = 2010K and 
£c/£o = 1-40. They measure a mean mass loss rate of 
10~ 4 Mq yr _1 during more than nine periods, each 400 
years long in their model. In one such period the star (in 
this case IRC +10216) would loose 0.04 M 0 and during 
the whole computation more than 0.36 M 0 , which in turn 
is more than half of the stellar mass M*. Therefore, the 
modelling of these long-term variations require a perma¬ 
nent “refilling” of the mass in the computational domain 
to prevent the model from changing due to the depletion 
of gas by the stellar wind. At present, our models do not 
permit such a “refilling” of matter, e.g. by transport across 
the inner boundary, as done by SID01. Therefore the cal¬ 
culations have to be stopped before a significant fraction 
of the mass contained in the computational domain is lost 
(typically of a few 10 2 years) and no studies of long-term 
variations can be performed at the moment, excluding a 
direct comparison with the result of SID01. Turning to 
short-term variations, since we find the stellar parameters 
above to be unrealistic in the sense of the very low value 
of the effective temperature we have not made a model 
with these parameters. 

Next, we examine the response of the model to various 
forms of the momentum transfer term. 
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Table 3. Detailed study of the momentum transfer term. 
All models below are drift models and have the same stel¬ 
lar parameters as model C19dl in Table [|. The fraction of 
specular dust-gas collisions are defined by e (see Sect. 3.2). 


model 

Cd (approx.) 

e 

(M) 

[M 0 yr x ] 

(Uoo ) 

[kms -1 ] 

C19dl 

s^tLA 

1.0 

1.4 10 -5 

24 

CHV 

smHV 

1.0 

1.110“ 5 

22 

CdO.O 

/''/LA 

0.0 

1.3 10“ 5 

23 

Cd0.5 

✓''/LA 

0.5 

1.4 10“ 5 

23 


4.3. Detailed study of the momentum transfer 


As mentioned in Sect. 3.2 several different forms of 


the momentum transfer term are in use in existing 
stellar wind models. In all previously discussed calcu¬ 
lations we have used the limits approximation of the 
drag coefficient, Cjj A , assuming fully specular collisions 
(e = 1.0). In this section we study the behavior of the 
C19dl model when we on the one hand use an other 
approximation of the drag coefficient, and on the other 
hand carry out the calculations using diffusive gas-dust 
collisions. Table || gives the different model configurations. 

The momentum transfer term makes the RHD3-system 
we describe very stiff. If the system is solved explicitly 
there might be a reason to use an analytically simple ex¬ 
pression for the drag coefficient such as Cp V (Eq. p4| ). 
Since this approximation is in use we find it useful to 
compare model C19dl, which is calculated using C^ A , 
with model CHV calculated using Cq V . This is specifi¬ 
cally motivated if we consider the measured limits of the 
speed ratio which we found for the models in the previous 
section (Eq.|25|). The two models are shown in Fig. 

By definition the drag coefficient Cq V is not suitable 
in regions of low drift velocities where its use introduces a 
systematical error in the solution. This observation is con¬ 
firmed in the diagrams of the speed ratio and the relative 
error of the drag coefficient (i.e. the drag force) used in the 
calculations, Fig. ||c,d. The speed ratio Sb in model CHV 
is systematically greater in the shocked regions where that 
of C19dl is low. As expected the drift velocity is the low¬ 
est in the regions behind dense shocks where the radiation 
pressure on the other hand is the largest. As discussed be¬ 
fore it is in these regions where a major fraction of the 
dust is located. 

Despite this it is not evident from the figure, and the 
values in Table that the global characteristics differ be¬ 
tween CHV and C19dl. We want to point out that we 
cannot be certain as whether a better numerical advection 
or other stellar parameters will change this conclusion. In 
view of these results and the fact that the speed ratio does 
not reach high numbers we conclude that the high-velocity 
approximation of the drag coefficient should be used with 
caution in AGB wind models showing similar speed ratios. 



10 


15 

R* 


20 


25 


Fig. 3. Study of the high-velocity drag coefficient approx¬ 
imation. The models are C19dl (solid line), and CHV 
(dotted line): a gas velocity; b gas density; c speed ra¬ 
tio; d relative deviation of the drag coefficient used in the 
calculations as compared to Cb; e relative error of the 
structure from being CMC, the drag coefficient used in 
the drag force is the expression used in the calculations, 
i.e. for CDHV and C^ A for C19dl. The drag force 
is volume integrated. The features in the three lowermost 
panels of model C19dl associated with the low values in 
the speed ratio at about 9.50, 19.0 and 23.5 i?* should not 
be taken into account in physical consideration (cf. Fig. |j). 
The same argument concerns model CHV. Note that the 
relative error of Cq V is significant in most of the envelope, 
especially in regions of lower speed ratios. 
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Fig. 4. Model response to two different types of collisions. 
100% diffusive collisions are assumed in model CdO.O (solid 
line). In model C19dl (dash-dot-dotted line) they are in¬ 
stead 100% specular: a dust density; b dust to gas tem¬ 
perature ratio; c ratio of the diffusive term Cb.dm to the 
common term Cb,common in Cb (Eq. [II]); d drift velocity. 
The diffusive term Cb.diff makes up a large part of Cb 
(up to ~ 40%). Still the overall behavior of model CdO.O 
largely agrees with model C19dl. 


Figure ^3 shows the relative difference of the drag force 
and the radiative pressure. This quantity can be used 
as a qualitative measure on how close CMC is achieved 
throughout the envelope. C19dl is close to CMC every¬ 
where (< 20%) except in the shock fronts at about 9.50, 
19.0 and 23.51?* where the dust velocity is lower than 
it would be in CMC. The same argument is valid with 
the shock fronts in model CHV at about 13.0, 16.5 and 
20.51?*. 

We now compare models with different values on e. 
We find that the difference in the drag coefficient between 


an expression describing fully diffusive and fully specu¬ 
lar collisions under the current conditions is small. The 
additional term used with diffusive collisions (Cb.diff) is 
proportional to the fraction of such collisions. As in the 
case discussed above, this term becomes important at low 
drift velocities in the shocks. In Fig. [] the fully diffusive 
collisions model CdO.O (solid line) is compared with the 
fully specular C19dl (dash-dotted). 

The diffusive term contains the temperature ratio of 
the dust and the gas temperatures, Fig. []b. The dust tem¬ 
perature is equal to the radiative temperature in AGB 
wind models using gray opacity and assuming radiative 
equilibrium. Likewise the gas temperature is mostly close 
to the radiative temperature. In Fig. Qc we see that the 
diffusive term (CD.diff in Eq. ^oj) is up to 40% in size com¬ 
pared to the common term (Cb,common)- Again the dif¬ 
ference is the largest in the regions behind shocks. The 
dust density and drift velocity plots (Fig. |]a,,d) show that 
the apparent effects on the structure in CdO.O compared 
to C19dl are minute. The dust likely needs to be more 
accurately handled by increasing the numerical precision 
in the advection calculation to second order and using 
drift related effects in the dust formation (KS97) before 
we can trace differences in dust quantities caused by the 
type of collisional distribution. We finally note that the 
differences will be even smaller with a lower fraction of 
diffusive collisions. 

5. Conclusions 

In this paper we have studied details of the dynamical in¬ 
teraction between dust and gas in cool stellar winds. This 
interaction is vital for the formation of a dust-driven out¬ 
flow. In our work we have studied changes of general prop¬ 
erties in the wind with the additional freedom of a separate 
equation of motion for the dust component, and we have 
not specifically attempted to model periodic winds. 

We have here first carried out a detailed study of the 
physical and numerical issues associated with a separate 
equation of motion including mass, momentum, and en¬ 
ergy interaction terms. The results of the drift models have 
been compared with position coupled (PC) models (i.e. 
models without the separate dust equation of motion). 

The results of the comparisons have been divided into 
those of local character, i.e. modified spatial distribution 
of the dust component, and global such as time-averaged 
mass loss rates and the general wind behavior. Of the local 
results we have found that the dust is more concentrated 
to shocked regions in drift models; the inter-shock regions 
are more depleted of dust than in PC models. Globally, we 
have found that the mean mass loss rates of the drift mod¬ 
els are about the same when compared with PC models 
for the massive winds. We have however not been able to 
generate the weakest winds, i.e. those that were station¬ 
ary in the PC models. The collisional gas-dust coupling in 
the atmosphere of these model configurations is not strong 
enough to form a dust-driven stellar wind, and the dust 
formed instead leaves the star alone. 
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We have found that a high-velocity approximation of 
the drag coefficient in the momentum transfer term (used 
in wind models by e.g. SID01) shows a systematical devi¬ 
ation in the drag force (compared to when using the full 
expression) throughout most of the envelope. The high- 
velocity approximation should be used with caution in 
winds such as those that we have modeled where the drift 
speed (and the speed ratio) is low in the regions where 
the most of the dust resides. One of the other two ap¬ 
proximations presented, Eq. (|23|) or Ecp (p2|), should be 
used instead. The winds of drift models calculated with 
diffusive collisions have turned out to be very similar to 
winds calculated with specular collisions. The interpreta¬ 
tion should, however, not be that diffusive collisions do not 
make a difference. This result may turn out to be very dif¬ 
ferent if frequency dependent opacities are used (instead 
of gra y), whic h modify the temperature structure in the 
wind ( Hofner 199S| ; [Hofner et al. 2002a| , bj) , and in partic¬ 
ular the dust/gas temperature ratio (see Eq. (To) ). 

The three-component wind model we have presented 
leaves some issues open that need to be addressed be¬ 
fore we will attempt to match observed stellar wind prop¬ 
erties with drift models. The first issue is numerical. In 
order to handle the calculation of the advection more ac¬ 
curately, a second order advection scheme instead of the 
current first order is necessary; our findings on numerical 
convergence properties discussed in this paper (relate d to 
the two-phase shock interaction mentioned in Sect. 3.2) 
could possibly also be moderated this way. Furthermore, 
we have so far neglected the fact that the dust formation 
is dependent on the drift velocity between the condensing 
material and the dust grains (e.g. KS97). An improvement 
in the description of the formation processes is thus phys¬ 
ically justified. A gas-dust interaction term we have not 
included in the current study is the heating due to fric¬ 
tion (see Sect. 2.3.2). We plan to discuss the importance 
of this term in more detail by including it in the calcula¬ 
tions. Most AGB stars are long-period variables, and the 
effects of stellar pulsations should be included. This will 
be done in the next generation of wind models. 
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Appendix A: Physical conditions in the wind 

In this subsection we comment on the justification in us¬ 
ing the continuum approximation for the gas component 
in the extended atmosphere of the AGB star. We also show 
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Fig. A.l. Knudsen numbers of the PC model C19c2. The 
presented instant shows an evolved stage of the wind: a gas 
Knudsen number, /C„; b gas-dust Knudsen number, fC' n d ; 
c gas density; The gas is safely in the continuum regime 
as /C n <C 1. The simplified gas-dust Knudsen number IC n d 
shows that the dust is in the free-molecular regime com¬ 
pared to the gas. These are both essential requirements 
in the derivation of the expression of the drag force used 
here. 


that the dust is in the free molecular regime when consid¬ 
ering gas-dust collisions. Both are essential requirements 


in the derivation of the drag force in Sect. 2.3 


Both the gas and the dust are dilute in the extended 
atmopsheres of the AGB star. Thus, 


> d 


(A.l) 


where S is the mean molecular distance; n is the num¬ 
ber density of the respective particles; d is the diameter 
of the gas or dust particles respectively. The continuum 
approximation is valid for small Knudsen numbers (/C n ) 
that satisfy, 


/C n = A gg /L -C 1 


(A.2) 


where A gg is the mean free path of a gas particle; L is the 
characteristic size of the body or system to which the flow 
properties are related. A scale length of a macroscopic gra¬ 
dient, such as that of the density, yields a better estimate 
on L than a constant (B94); 



(A.3) 


An order of magnitude estimate of the mean free path, 
(assuming an elastic collision cross section er g ) is given by 

Agg = { nf7 g) ■ (A-4) 


Typically the collisional cross section in a gas is about 
<T g = 10“ 15 cm 2 . The Knudsen number as defined in 
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Eq. (A.2) is shown in Fig. |A.l| a for the PC model C19c2. 
The Knudsen number clearly stays well below 1 through 
all the extended atmosphere we model, hence we conclude 
that the use of the continuum approximation is justified. 

The size of the dust grain (d d ) is the body charac¬ 
teristic length L in dust grain-gas particle collisions. The 
mean-free path Ag d of a gas particle that has (just) collided 
with a dust grain, that drifts through the gas (So > 0 ), 
may be much shorter than that of the rest of the gas par¬ 
ticles (B94). An upper limit of the associated Knudsen 
number assuming zero drift speed (So = 0 ) is given by, 


^n,d — ^gg/^d 


(A.5) 


IC n d is shown in Fig. ATb, where we see that it is > 10 10 


all through the atmosphere. The mean free path A( d un¬ 
likely drops by a factor of ten when the speed ratio stays 
within the limits 0.1 < So < 10 (see Eq. |25|). We con¬ 
clude that the dust grain is in the free molecular flow 
regime compared to the gas. In the current context a free 
molecular flow means that the velocity distribution of a 
gas particle incident on a dust particle is independent of 
the history of previous collisions (involving other gas par¬ 
ticles) with the same dust particle. This last statement is 
however not equivalent as to consider the dust by itself 
to be in the free molecular regime. Such a conclusion re¬ 
quires a measure of the dust-dust collision mean free path; 
a measure which we have not made in this paper. 
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